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SUMMARY 


In a recently published work by Abarbanel and Gottlieb (NASA CR-1 59386 ) , 
a new class of explicit time-split algorithms designed for application to the 
Navier -Stokes equations for compressible flow was developed- These algorithms, 
which utilize locally-one-dimensional (LOD) spatial steps, were shown to possess 
stability characteristics superior to those of other time-split schemes- In 
the present work, the properties of an implicit LOD method, analogous to the 
Abarbanel-Gottlieb algorithm, are examined using the two-dimensional heat con- 
duction equation as the test problem. Both temporal and spatial inconsistencies 
inherent in the scheme are identified. The principal result of the present work 
is the development of a new consistent, implicit splitting approach. The 
relationship between this new method and other time-split implicit schemes is 
explained, and stability problems encountered with the method in three dimen- 
sions are discussed. 


INTRODUCTION 

Many of the methods currently in use for numerically solving the Navier- 
Stokes equations for compressible flow are based on the concept of time-split 
and fractional-step finite-difference schemes developed at length in reference 1 . 
These schemes include the Douglas-Gunn type alternating-direction implicit (ADI) 
methods of Briley and McDonald (ref. 2) and Beam and Warming (ref. 3) and the 
time-split explicit approaches such as that of MacCormack and Baldwin (ref. 4) . 
The primary impetus in the development of time-split schemes is to reduce the 
amount of computational work to advance the solution one time step. In an 
implicit method, direct inversion of the full matrix required to gain the solu- 
tion of the algebraic system associated with the difference equations requires a 
prohibitive amount of computational effort. By (approximately) splitting the 
matrix, this computational effort can be substantially reduced. In explicit 
methods, the splitting of the equations so that the various operators are 
advanced separately allows one, in principle, to advance each of these steps 
at its own stability limit. Since the stability limit for some of the steps 
can be substantially larger than for others in a typical high Reynolds number 
flow, a savings in computational effort can be realized. Currently, available 
splitting schemes suffer from various shortcomings. The spatially split 
Douglas-Gunn type methods (e.g., refs. 2 and 3) have stability problems in 
three dimensions, while the time-split explicit method of reference 4 does not 
achieve the sought-after independence of stability criteria for the various 
split steps. 

Recent work by Abarbanel and Gottlieb (ref. 5) has shed new light on the 
stability restriction of the time-split scheme of MacCormack and Baldwin 
(ref. 4) and opened up the possibility of a new class of implicit methods. In 
the current paper, the consistency of time-split schemes as introduced by 
Abarbanel and Gottlieb is examined, and a new consistent time-split finite- 
difference algorithm is introduced. 



In the time-split scheme of reference 4, the mixed-derivative term for 
viscous flow in the Navier-Stokes equations is apportioned among the various 
space-split operators. Abarbanel and Gottlieb point out, through a stability 
analysis of the full (linearized) Navier-Stokes equations, that this does not 
constitute an optimal split in that the allowable time step for each split step 
is influenced by the mesh spacing associated with other split steps. In ref- 
erence 5, Abarbanel and Gottlieb propose a new method of splitting in which the 
spatial operators of the Navier-Stokes equations are first split into hyperbolic 
(i.e., Euler equations), parabolic, and mixed-derivative operators, and then 
each of these operators is solved with a time-split explicit scheme. The method 
is thus akin to the locally-one-dimensional (LOD) schemes, i.e., schemes in 
which only one term in the spatial operator appears in each split step. This 
scheme is proven to be an optimally split scheme. Further work by Abarbanel 
and Gottlieb indicates that such an operator-split scheme, in which the LOD 
explicit-difference method is replaced by an LOD backward Euler implicit method 
(with the cross-derivative operator advanced explicitly) , is unconditionally 
stable and possesses a very good smoothing rate (a desirable feature for the 
multigrid method) . Such a method thus appears to offer an attractive alterna- 
tive to the spatially split Douglas-Gunn schemes of references 2 and 3, which 
are known to have stability problems in three dimensions (ref. 6) . 

As a summary of the status of current finite-difference methods for solu- 
tion of the three-dimensional Navier-Stokes equations, the desirable properties 
of an "ideal" scheme are listed; 

1 . Unconditional stability 

2. Consistency 

3. Temporal accuracy 

4. Ability to recover steady-state solutions, when they exist, independently 

of the iteration history 

5. Ability to perform required matrix inversions at a minimum of computa- 

tional time 

The currently used Douglas-Gunn methods of references 2 and 3 meet criteria 2 
through 5. On the other hand, the implicit version of the scheme developed by 
Abarbanel and Gottlieb in reference 5 has been demonstrated to meet criterion 1 , 
but its properties in regard to criteria 2 through 4 are unknown. 

The motivation of the current investigation was to understand better the 
properties of the LOD-type algorithm as proposed in reference 5. The test prob- 
lem studied was the solution of the heat conduction equation on the unit square 
subject to steady Dirichlet data. Numerical solutions to this problem using both 
explicit and implicit LOD schemes revealed two immediate difficulties: (1) a 

temporal inconsistency in the steady-state solution due to the operator split- 
ting and (2) errors in the solution due to the application of boundary data. 

The first difficulty is well known and there are a variety of ways of allevi- 
ating it. The second difficulty is not so well known and is somewhat more sub- 
tle but has been studied to some extent in the Russian literature (ref. 1 ) and 
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is discussed briefly by Mitchell (ref. 7) . The errors produced by the temporal 
inconsistency are of the order of the temporal truncation error and, hence, pro- 
vide negligible contamination of the solution with explicit schemes, since sta- 
bility considerations dictate that this time step be bounded below the spatial 
truncation error. However, stability requires no time-step bound for implicit 
schemes, so that the temporal inconsistency can actually swamp the true solution 
for large values of the time step. Computational verification of these errors 
for both explicit and implicit LOD schemes is presented in this paper. 

As discussed in references 1 , 7, and 8, the boundary-condition errors in 
the LOD methods arise frcxn the fact that the intermediate solutions do not 
represent consistent approximations to the dependent variables at any time 
level. Hence, the intermediate boundary data may not properly be the given 
boundary data for the problem and, in general, will be dependent on the form 
of the split operator and the given boundary data. Consistently split schemes 
do not suffer frcan this problem, since all intermediate solutions represent 
approximations to the dependent variables at the new time level (ref. 8) . The 
method of undetermined functions (MUD) developed in reference 1 is shown to 
alleviate the boundary difficulty for the LOD schemes for Dirichlet data. A 
generalization of this procedure leads to the construction of a general consis- 
tent LOD scheme which requires no boundary-error correction for steady data. 

This new approach is presented in the current work. 

In the remainder of this paper the inconsistencies associated with LOD 
schemes are developed in detail and methods for their alleviation outlined. 

The information thus developed is used to devise a general consistent LOD 
scheme. Comments are then made on the relation of this scheme to other ADI 
schemes. 


S!iMBOLS 

A£ full-step tridiagonal-difference operator defined in equation (1 3) , 

fi. = 1 ,2 

Aj^ half-step tridiagonal-difference operator defined in equation (30) , 

= 1 ,2 

% partial-step tridiagonal-difference operator defined in equation (62), 

A = 1 ,2 

Bj, full-step explicit-difference operator defined in equation (1 4) , 

I = ^ ,2 

dji boundary-condition residual function for LOD forward Euler scheme 

defined in equations (A3), A = 1,...4 

D discrete computational domain (fig. ^) 

D^ continuum domain (fig. 1) 
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D-| ,I>2 


union of all horizontal and vertical grid lines in respectively 
(fig. 1) 

D] ,D 2 union of all horizontal and vertical grid lines in D u T, respec- 
tively (fig. 1 ) 

boundary-condition residual function for LOD Lax-Wendroff scheme 
defined in equations (A7) , il = 1 , • . • 4 

E error in solution defined as difference between analytical and 

numerical solution 

f£ boundary-condition function, 9. = 

F£ constants in exact solution to test problem, = 1 ,2 

g£ boundary-condition residual function for LOD backward Euler scheme 

defined in equations (33) to (36), A = 1,...4 

gj^ boundary-condition residual function for LOD predictor-corrector 

scheme, 9^ = 1 , ... 4 

g£ boundary-condition residual function for LODQ scheme defined in 

equations (69), 9. = 1,...4 

G Von Neumann damping ratio 

h spatial step size (subscripts indicate spatial direction) 

Hj horizontal grid line in D, j = 1,2,...J 

i = ^ 

I total number of spatial grid points in x-| -direction 

J total number of spatial grid points in X 2 -direction 

k integer exponent 

L-| ,L 2 /L 3 general time-split spatial-difference operators 
M*| ,M 2 fM 3 general LOD split-difference operators 
n time-increment counter 

q correction function for LOD scheme determined by method of 

undetermined functions 

Q correction function for LODQ scheme 

R steady-state residual defined in equation (1 7) 
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R 


steady-state residual of whole step form of time-split schemes 



u dependent variable 

Vi vertical grid line in D, i = 1,2, ...I 

x*|,X 2 ,X 3 Cartesian coordinate directions 
a parameter that sets temporal order of accuracy 

= 2 At h“2(l - cos 0£) , = 1 , 2 

element of discrete boundary F, il = 1 , . . . 4 
element of continuum boundary F^, il = 1 , • . • 4 
F discrete boundary (fig. 1) 

F^ continuum boundary (fig. 1 ) 

6u/6t general time-difference approximation 

6u/6x^ spatial-difference approximation, il = 1,2,3 
At time step 

ri£ = At h“^ sin 0£, il = 1 , 2 

0j^ phase angle of Fourier error analysis, il = 1 , 2 

Aj^ second-order central second-difference operator defined in equa- 

tions (6) and (7) , il = 1 , 2 

T normalized time step defined in equation (20) 

element of the set il = 1 , • . • 4 

SI set of grid lines one mesh increment interior to discrete boundary 

F (fig. 1) 

Superscripts: 


n 

time- step 

counter 

m 

time-step 

counter 
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(m) 


time-split-step counter 


* intermediate time-level indicator 

Subscripts: 

i x-| -mesh increment counter 

j X 2 -mesh increment counter 

max maximum 

min minimum 

1,2 Cartesian coordinate direction indicator 

Abbr e V i at i ons : 

ADI alternating direction implicit 

DG Douglas-Gunn 

LOD locally-one-dimensional 

LODBE LOD backward Euler 

LODBEC consistent LOD backward Euler 

LODFE LOD forward Euler 

LODPC LOD predictor-corrector 

LODQ LOD with correction function Q (eqs. (60) to (66)) 

MUD method of undetermined functions 

MUDB MUD with boundary data adjusted 

MUDF MUD with correction function applied over entire domain D 

2- D two-dimensional 

3- D three-dimensional 

DEVELOPMENT AND APPLICATION OF SIMPLE LOD SCHEMES 
MODEL PROBLEM FORMULATION 

The model problem chosen for study in this paper is the heat conduction 
equation in two dimensions on the unit square. That is. 
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( 1 ) 


9u 

3t 


02u 92 u 

a 2 ^ o 2' 

3x, axj 


will be solved 
Dc = 


on the domain which is the open set 

^{x-j ,X 2 ) : 0 < x-| ,X 2 < l} 


The boundary of the continuum domain is the union of the four unit-length 

line segments YcJl^ illustrated in figure 1 . The discrete domain D on 
which the computations are carried out is obtained by overlaying with a 

uniform net, with D defined by 



where I and J are the number of mesh points in the x*| - and X2-directions, 
respectively, and the mesh spacings are given by 


1 



1 



In all numerical solutions performed in this paper we will take I = J. Conse- 
quently, the spatial step sizes are equal (i.e., h*| = h2) # so that the unsub- 

scripted symbol h is used to denote both spatial increments. The discrete 
boundary of D, the set F, is the union of the four discrete line segment 
sets Yilf which have obvious definitions. The discrete region is depicted in 
figure 1 . To augment the results which follow, the discrete set ^ is defined 
as the union of the four line segments lying one mesh width frcxn F within 

D (fig. 1) - 

As a final point, additional notation is introduced to simplify the 
description of the manner in which the algorithms described in this paper are 
performed computationally. Let 
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These two sets are illustrated in figure 1 . Note that Hjj, is the S/th 
horizontal line and VJ, is the itth vertical line in D ^ F. Now/ define 


J -1 

Di = U HJl 
A =2 


J 

D, S U % 
1=1 


1-1 

D2 = U Vjt 

1=2 


I 

D2 = U VS, 

il =1 


The sets D] and D2 are equivalent to D but merely viewed in another way - 
as the union of all horizontal and vertical lines in D, respectively. The sets 
D-j and D2 include the adjacent boundary lines. 

During the course of this paper, the numerical solution to equation ( 1 ) 
using various algorithms is discussed for two separate problems, designated test 
problem ^ and test problem 2 (or TPl and TP 2 ) . Both problems have Dirichlet- 
type boundary data. Let 


u(XT,X2,t) = fs,(XT,X2,t) (X] ,X2) ^ Ycg,f ^ 


(2) 


The boundary data for the test problems are 


f-| (x-| ,X2/t) = F-] sin HX2 

f2(x-],X2ft) = F2 sin nx] 

f3(xi,X2ft) = F-[ sin Ux2 

f4(x-],X2ft) = F2 sin TTxi 


J 


( 3 ) 
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where F-| and F 2 may take on the values 0 or 1 . The steady-state analytic 
solution to equation (1 ) on subject to the boundary data given by equa- 

tions (3) is 


u(xi,X2) = F2 sech cosh n^X2 - 


sin TTxi 




cosh 

"(«i - 2) 


sin 7 Tx2 


(4) 


For TPl , F-| = 0 and F 2 = 1 ; for TP2, F-| = 1 and F 2 = 0. The boundary con- 

ditions and applicable steady-state analytic solutions for TPl and TP2 are plot- 
ted for illustrative purposes in figure 2, 

The reason for the introduction of the two separate test problems arises 
fran the nature of the errors in time-split difference methods illustrated in 
this report. As discussed herein, errors are introduced in time-split methods 
as a result of inconsistencies associated with the boundary conditions. The 
particular boundaries associated with these errors vary depending on the par- 
ticular time-split scheme, hence the need for the two problems. 


DIFFERENCE APPROXIMATIONS 

All numerical schemes discussed in this paper use first-order, one-sided 
(backward or forward) finite differences for the time derivatives and second- 
order, central differences for spatial derivatives appearing in equation (1 ) . 

If U£j is the approximation to u(ih'|,jh2rn At), then the difference forms 
are 


\n-D ^ 

/ ^ ^ij “ ^ij 
\ 3 t/ij At 


(5) 


where 0 S p ^ 1 and At is the time step. The spatial derivatives are 
approximated by 


Ux2y.. h2 

\ 1/ij 1 


( 6 a) 


= A, 



( 6 b) 
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H, j+1 


- 2u 


ID 


■ n 


^2 


(7a) 



(7b) 


When no confusion can arise regarding a given spatial location, the notation 
given in equations (5) to (7) is simplified by dropping the subscripts on u. 
For example, the forward-time, centered-space approximation to equation (1) 
may be written as 


un - un-1 

= (At + A2)un-1 


( 8 ) 


with the implied understanding that this approximation is to be applied at all 
points in D. 


TIME-SPLIT SCHEMES 
General Idea 

Time-split schemes are methods in which the computations associated with 
advancing a numerical approximation a time distance At are divided into a 
sequence of two or more substeps. This division is done to obtain an algorithm 
with less ccxnputational work than the whole-step scheme. Time-split methods 
may be implicit, explicit, or a combination - sane implicit steps, sane explicit 
steps. Canbination-type schemes are often referred to as predictor-corrector 
or hybrid methods. 

As an illustration of the splitting concept, consider a general two-step 
scheme for equation (1 ) which advances time from level n to n + 7 . Such a 
scheme may be written in operational form as 



L-j (A-| u,A2U) 


(9a) 



L2(A-|U,A2u) 


(9b) 
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where — - represents sOTie approximation to — , and Li and Lo are specified i 
6t 9t 

functions which define the nature of the splitting. The variable u with no 
time-level indicator is meant to imply that a number of time levels may be 
associated with each A^,. Combining equations (9a) and (9b) gives the whole- 
step scheme, which is equivalent to the split scheme at interior mesh nodes 
(i.e., for points in D) . The whole-step algorithms assume the general form 


un +1 _ yn 

At 


(A) 



a)u^ + 


+ At*^L3(A] u,A2U) 


(9c) 


where k is an integer, and L 3 is an error term resulting from the split 
and depends explicitly on L-| and L 2 - Equation (9c) may not be valid on or 
near the boundary of D (e.g,, on one of the lines in Q) . For locally-one- 
dimensional (LOD) splittings, each of the split operators (eqs. (9a) and (9b)) 
involves spatial derivatives in only one direction. Therefore, the general form 
of the time-split LOD schemes for equation (1) may be written as 



1 


Ml (Aiu) 


( 10 a) 



M2 (A2U) 


( 10 b) 


with the corresponding whole step 


un+1 _ ^n 

At 



- a)u^ 


+ au 



+ At*^M3 (Ai u,A2u) 


( 10 c) 


The Douglas-Gunn Method 

One of the most widely used implicit time-split methods is that due to 
Douglas and Gunn (ref. 9). For the scalar equation (1), the spatially split 
Douglas-Gunn (DG) scheme is 


u 


* 


At 


= Aiu* + A2U*^ 


( 11 a) 
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( 11 b) 


un+1 _ uii 

It 


= A-j u* + A2U*^‘*'^ 


The asterisk is a symbol of convenience and can denote some intermediate time 
level/ say n At < t* ^ (n + 1 ) At. As shown in this section, however, for some 
time-split methods, u* or other intermediate levels may not be a consistent 
approximation to the solution of equation (1 ) . Note that for the DG scheme each 
split step (eqs. ( 11 )) is a consistent approximation to equation ( 1 ). 

From an operational standpoint, equations (11) may be rewritten as 


A-| u* = B 2 U^ 

(XI ,X 2 ) 

A 2 U*^"'‘^ = u*^ + At A-| u* 

(XI ,X 2 ) 


where 


AS, = (1 

- At Aj^) 

II 

Bs, = (1 

+ At Aj,) 

il = 1 ,2 


e Dt 

( 12 a) 

e d 2 

(1 2 b) 


(13) 


(1 4) 


The equivalent whole-step DG algorithm is obtained as follows. Multiplying 
equation (1 2 b) by A*; yields 


AtA 2 u'^''’^ = At + At A] (A 7 U*) (15a) 

= A-iu'^ + At A-| (B 2 U*’) (15b) 

Note we have commuted A-| and A-] to obtain equation (1 5a) and used equa- 
tion (12a) to obtain equation (15b). If equation (15b) is expanded using the 
definitions given by equations (13) and (14), the equivalent DG whole-step 
fonnula in difference form is 


un+1 _ un 

At 


(Al 


+ A2 )u«+^ 


- At^A-| A2 




(XT ,X 2 ) e D 


( 16 ) 


The time-splitting error term for the DG method is formally second order in time 
and, in addition, is proportional to an approximation to the time derivative 



Uf If the problem being solved has a steady-state solution/ the time-splitting 
error term would be expected to vanish as the steady solution is approached. \\ 

However, this will not be true of an unsteady problem. The symbol R is used 
to denote the so-called steady-state residual: 


R = (At + A2)u 


( 17 ) 


The entire residual for this and other schemes is denoted R. For the DG method, 
R is given by 


Rn+l 


(A-j + - At2A-|A2 



0 8) 


The DG method is unconditionally stable, as may be verified by a simple 
Von Neumann analysis. The damping ratio G is given by 


1 + 3 ] 32 

(1 + 3 i ) (1 + 32 ) 


( 19 ) 


where 3^ = 2 At h ^(1 - cos 0£) • 

The DG approach was applied to both TPl and TP2. The method is implemented 
by first implicitly solving equation (1 2a) (via the Thomas algorithm) along all 
interior horizontal grid lines in D (i.e., D-| ; see fig. 1). This is followed 

by a similar implicit solution of equation (1 2b) along all vertical lines in D 
(i.e., D 2 )« The time splitting thus has the computational effect of approximat- 

ing the inversion of the original block tridiagonal matrix (frcan (A-j + A 2 )u^"‘“^) 
by a sequence of simpler scalar tridiagonal inversions of the unidirectional 

operators. The penalty for use of such an artifice is, of course, the time- 
splitting error. The computational results are summarized in table I for three 
different time steps. The normalized time step T is defined by 


T = 



( 20 ) 


where t = 1 corresponds to the stability limit of the unsplit forward-time, 
center ed-space method of equation (8) . Table I gives the maximum and minimum 
absolute errors IeI^^^x 1^1 min maximum steady-state residuals R 

were driven to machine zero, approximately 10“1^), which demonstrate the 
h^-convergence of the computed solution to the exact solution. In these calcu- 
lations the maximum errors occur at the center of the grid, as should be 
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expected since this location is farthest frcxn the boundaries where the solution 
has no error. 


Locally-One-Dimensional Methods 

In this section three locally-one-dimensional methods are discussed: 

(1 ) the LOD forward Euler method, (2) the LOD backward Euler method, and (3) the 
LOD predictor-corrector method. 

LOD forward Euler scheme .- The LOD forward Euler (LODFE) scheme is a two- 
step time-explicit approximation to equation (1 ) . The two steps are 


u* - u^ 

= A*] u^ 

At ' 


(21a) 


un+l _ u* 

At 


= A2U* 


(21b) 


Using equation (1 4) , the computational steps are 


I* = B-]U^ 


(XI ,X2) s D] 


u* = f£ 


(X-] ,X2) ^ = 1 , 3 


(22a) 


^n+l = B 2 U* 

(X] ,X2) 

e D2 

n -t-1 i-ri"*"1 

(jH+l = 

(XI ,X2) 



(22b) 


* n+1 

For TPl and TP2 we take f£ = fj^ = fj^, for (x-| ,X 2 ) ^ Yjlf ^ = 1,...4. It is 
easily verified that the stability restriction is T ^ 2. The equivalent whole- 
step method is obtained by substituting equations (22a) into equations (22b) 
and expanding the operators: 


jjn+1 « ^n 

= (An + Ao)u^ + At A-iAou^ (x*|,X2) ^ D (23) 

At 

Note that the time-splitting error term At A*|A 2 U^, unlike that for the DG 
scheme (eq. (16)), is first order in At and depends upon the function itself 
u^ rather than a form of u^. 


The solution, error, and residual maps generated by the LODFE algorithm 
applied to TPl are presented in figure 3 for t = 1 and x = 2. Note that the 
maximum and minimum errors now depend upon At, as predicted by equation (23) , 
and that the steady-state residual (eq. (17)) is no longer arbitrarily small 
as with the DG scheme. The error map in figure 3 introduces another interesting 
point - the maximum error and residual now appear along grid lines immediately 
adjacent to a boundary, that is, along elanents of the set (See fig. 1.) 

LOP backward Euler scheme .- The two-step implicit LOP backward Euler 
(LOPBE) method for equation (1 ) is given by 


u* - u'^ 

At 


= Ai u* 


(24a) 


un+1 _ u* 

aT" 


= A 2 U 


n+l 


(24b) 


with ccanparable computational steps 


A-) u* = u*^ 


u* = fa 


(x-| ,X 2 ) s Pi 


(XI ,X2) e Yi,, = 1,3 


(25a) 


A 2 U 


n+l = u* 


(xi ,X2) s P2 


uH+1 = f?+^ 


(XI ,X2) e Y)l» *• = 2,4 


(25b) 


For TPl and TP2 we take fj^ = for (x*| ,X2> ^ ^ = lr ...4 and 

is defined in equation (13). Multiplying equations (25b) by , substituting 
from equations (25a) / and expanding, we obtain the equivalent LODBE whole-step 
equation 


un+1 _ ^^n 

— = (Ai + A2)u^+1 - At A-|A2U^+^ (xi,X2) ^ D (26) 


Equation (26) indicates that the LODBE algorithm has a time-splitting error 
behavior similar to the LODFE scheme, namely, first-order temporal accuracy 
and a function-dependent splitting error. Computational results verify this. 


1 5 


I 



as shown by the error, residual, and steady-state solution maps for T = 1, 10, 
and 1000 presented in figure 4. Again, the steady-state solution is highly 
dependent upon T and the maximum error appears on the set adjacent to the 
boundaries. Note, however, that the implicit nature of the LODBE scheme has 
led to a desirable property: the method is unconditionally stable. 

LOP predictor-corrector scheme * Before proceeding to a detailed analysis 
of the errors that appear in LOD methods, we introduce one additional LOD 
algorithm, which Yanenko (ref. 1 ) calls the predictor-corrector splitting 
scheme and which we term LODPC. This approach combines a LODBE predictor with 
an explicit, leapfrog corrector in three steps as follows: 


u - u 


At/2 


n 


= A-| u^ 


(27a) 


u^*^^ /2 _ u* 

At/2 


A2U^’^^ /2 


(27b) 


un+1 _ ^n 

At 


(A-j + A2 )u^'*'^/^ 


(27c) 


The corresponding computational equations and whole-step scheme are 


A*| u* = 

(XI ,X2) 

e Dt 

* 



u* = f£ 

(X] ,X2) 

^ Yjl 

A2U^*^^/2 - u* 

(XI ,X2> 

e D 2 

un+l/2 = fO+V2 

(XT ,X2) 


uH"^^ = u^ + At (A-| + A 2 ) 

(XT ,X2) 

e D 

un+1/2 = 

(XI ,X2) 

^ YJ,. 


= 2,4 


\ (28a) 


(28b) 


\ (28c) 


= 1 ,...4 


J 


1 6 


ri n r 


in 


I III I 


II I 


III nil II 


ufl+l _ qH 


At 


= (A-| + A2) I 


'un+1 + u>A At 2 


4 At 




(XI ,X2) e D (29) 


★ n +1 /2 n +1 

For TPI and TP2 we take = fjj, ' = fjl ~ (xi/X 2 ) ^ Yjlf ^ = lr--*4, 

and 


^ ^ A. 


(30) 


The development of equation (29) requires the assumption of commutativity of 
the operator product A1A2 and (A-| + A2) - Note that the addition of the leap- 
frog corrector has had two effects on the whole-step formula as evidenced by 
equation (29) . First, the method is now second order in time; second, the 
splitting error now is proportional to u^- rather than u. The LODPC method 
is also unconditionally stable. The discussion of LODPC computational results 
is presented in the section on ccxnplete consistent schemes. 


SOURCES OF ERROR IN LOD SCHEMES 

For the sake of completeness, the sources of error in the simple forward 

and backward Euler LOD schemes are discussed here. A major source of the error 

in these methods arises from the well-known fact that the whole-step formulas 
(i.e., the scheme with u* eliminated) resulting from equations (21) and (24) 
are not consistent with equation (1 ) . Another error arises from a lesser known 
source, that is, from errors due to the application of the boundary data. These 
two sources of error are termed temporal and spatial inconsistencies, respec- 
tively, in this report. 

The spatial inconsistency for the backward Euler scheme (eqs. (24) ) is 

discussed in detail here following Yanenko (ref. 1). Examination of the struc- 

ture of this scheme reveals that calculations on the first sweep are influenced 
by boundary data applied along Yi and Y 3 but are in no way influenced by 
boundary data applied along Y 2 Y4* Similarly, on the second sweep the 

calculations are influenced by the boundary data applied along Y2 Y4 and 

not influenced by the boundary data applied along Yi and Y3- This fact is 
graphically demonstrated in figure 5. As shown, on a given sweep the effect of 
the spatial operator contained within the equation for the sweep renders the 
solution independent of data on any grid line except the one being swept. From 
equations (24), one also notes that as At ^ the coupling between the two 
sweeps vanishes. One is thus left with the solution of two independent series 
of boundary- value problems on the two sweeps. This fact is further reflected 
in the damping of the scheme given by 
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G = 


1 


(31) 


(1 + 3 i ) (1 + 32 ) 


where G is the ratio of the amplitudes of the Fourier modes between successive 
time steps. From equation (31)/ 


lira G = 0 


(32) 


which indicates that one approaches a direct solution to some (incorrect) steady 
problem as the time step increases. This behavior is in contrast to that of 
the Douglas-Gunn scheme/ in which the entire spatial operator appears in each 
sweep equation. (See eqs. (11).) Thus/ calculations on each sweep are influ- 
enced by data applied on all four boundaries. For the Douglas-Gunn method, 

G 1 as At 

The incOTipatibility of the two sweep equations (25) with the boundary data 
is illustrated by the fact that if equations (25a) is applied along the boundary 
segment Y2' 

^2 “ ^2 " A] f 2 = 92 (xi fX 2 ) ^ Y2 (33a) 

* 

Now, if f 2 is a function of time and we try to identify f 2 with some inter- 

n+1 /2 

mediate value of f, say f 2 r then equation (33a) becomes 


_n+l /2 _n 

f2 - f2 - 

At =92 

(x-| ,X 2 ) e Y2 

(33b) 

Note that, in general. 

92 does not vanish. 

For steady boundary data. 

92 is 

-At A-| f2 = 92 

(x^ ,X2) e Y2 


(33c) 


which again does not vanish, in general. Also note that 

92-0 (^1 ^ D u Yi u Y 3 u Y 4 (33d) 


1 8 


■ 


■ ■ II iiiiia « III 1 1 III! Ill im nil 


Similar relations along the other boundary segments may be written; 


£4 - £4 - At A-] £4 = g 4 
94 = 0 


{xi ,X2> € Y4 

(XT ,X2> e D.u Yl u Y2 u Y3 


( 34 ) 


- £* - At A2fi''^^ = 91 (xi,X2) e yi 


91 = 0 


(xi ,X2) s D u Y2 u Y3 u Y4 


(35) 


-n+1 * . . -ii+l . . _ 

^3 ■“ ^3 ■“ ^2^3 = 93 (^1 '^ 2 ) ^ Y3 


93 = 0 


) (36) 


(XI ,X 2 ) 6 D u Yl u Y2 u Y4 


With these definitions, equations (25) are now rewritten as 


A-| u* = u*' + (g 2 + 94 ) 


(xi ,X2) e Di 


A2u’^+'' = u* + (gi + g3) (xt ,X2) g T>2 


(37) 


If u* is eliminated from equations (37) , the whole-step scheme becomes 


A“|A 2 U^"*'^ = + (92 + 94 ) + A-| (91 + 93 ) (x-j ,X 2 > ^ D U T 


or, upon expandin 9 . 


un +1 _ un 

= (Ai + A2 )u^+' 

At 


At A-|A2U^**‘^ 


+ At-'' |j[g 2 + 94 ) + Ai (gi + ga)^ (xi ,X 2 ) e d u F (38) 

where A-|A 2 U^“*‘^ = 0 for {x-| ,X 2 > ^ F. Thus, the steady-state residual of the 
scheme 9 iven by equations (37) is 
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R = (A*| + A 2 ) u — At A 1 A 2 U 

+ At““^ [[(92 + 94) + A-| (91 + 93)j 


(XI, X2) e D u r 


(39) 


As noted in equations (33) to (36) , the functions 9 ^ are nonzero only on the 
boundaries. However, as illustrated in fi 9 ure 6, if 9 i 0 on Yi / then 

A*i 9i 5 ^ 0 on to-j 


At 91 = 


A] 93 = 


Thusr we would expect that the incompatibility of the split sweep operators with 
boundary data applied in the sweep direction contributes an error term to the 
steady-state residual of the whole-step scheme only on ^ (eq. (39) ) . This 
incompatibility contributes an error term to the residual R only on ^ 
because of the unique form of the functions 9 ^; however, the error in the solu- 
tion u caused by this term is propa 9 ated over the entire 9 rid. 

The errors present in the solutions in fi 9 ure 4 are a result of the addi- 
tional terms identified in the steady-state residual formula (eq. (39)); i.e.. 


and similarly on 003 . In fact. 


At 

“91 (xi .X2) ^ 0)1 


Elsewhere 


(40a) 


At 

— : 93 ^^ 2 ) ^ 0)3 

h2 


Elsewhere 


(40b) 


R = -At A 1 A 2 U + Af"^Ai ( 9 i + 93 ) (xi ,X 2 ) ^ D 


(41) 


as demonstrated in the next section. Note also from fi 9 ure 4 that the maximum 
error occurs on This arises frcxn the fact that the term At~^Ai (91 + 93) 

appearin 9 in equation (41 ) is nonzero only on 

Errors of the type discussed here are present in all time-split or LOD-type 
schemes. Similar boundary-error analyses are carried out for the LOD forward 
Euler, LOD Lax-Wendrof f , and a hybrid LOD scheme in appendix A. Procedures 
desiqned to eliminate these errors from time-split schemes are introduced in 
the next section of this paper. 
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COMPLETE CONSISTENT LOP SCHEMES 




REMOVAL OF TEMPORAL INCONSISTENCY 

A desirable property of difference schemes for the solution of equation ( 1 ) 
is that they be accurate for the unsteady case and recover any existing steady- 
state solution independently of At* Such a scheme is designated "complete 
consistent" in reference 1 • We now attempt to render the inconsistent LOD 
backward Euler method complete consistent through the removal of both the tem- 
poral and spatial inconsistencies. Removal of the temporal inconsistency is 
discussed first, since it is the most easily achieved. A similar modification 
of LODFE appears in appendix A. 

The LODBE method can be made temporally consistent by the explicit addi- 
tion of the term At^A-jA2U^ to the first step of the time-split operator. 
{At2A-|A2U*^ is defined to be zero for (x*| ,X2) ^ T.) Thus, the computational 
algorithm expressed by equations ( 37 ) becomes 

A-| u* = u^ + At^A*|A2U^ + (g2 + 94) (x-| ,X2) ^ D-j ( 42 a) 

;^2un+l = u* + (g-| + g3) (x*! ,X2) ^ D2 ( 42 b) 

where the boundary-condition residual functions have been retained for compati- 

bility with analyses which follow later. In whole steps the scheme (eqs. ( 42 )) 
is 


^n+1 „ ^n 
At 


(A-| + A2)u^'^'* - At2A-|A2 



+ Af"^ (gi + g3) + (g2 + 94)^ (xt ,X2) e d u T ( 43 ) 

Note the term that brings about the temporal inconsistency has been removed (for 
problems having a steady-state solution) by converting it to a u^ form. If 
we restrict our attention to D, g2 and g4 can be dropped from equation ( 43 ), 
since they are zero by definition on D, so that we obtain 


.n+l 


u' 


n 


At 


= (A-| + A2)u’^“*"^ - At^A*jA2 


^un+1 ^ ^n 

^ At > 


+ 




(x-| ,X2> ^ D 


( 44 ) 
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Recall that, in view of the analysis presented in the previous section, the term 
A*| (g^ + 93 ) cannot be eliminated, as it is nonzero at certain locations in D. 
(See eqs. (40)-) The method specified by equations (42) is termed the consistent 
LOD backward Euler method (LODBEC) . The technique used here to remove the tem- 
poral inconsistency was suggested to the authors by David Gottlieb, of Tel Aviv 
University. 

Another LOD scheme which is temporally consistent in its basic form - the 
LODPC method - was introduced earlier. If boundary-error problems are antici- 
pated, the computational scheme given in equations (28) may be rewritten as 


A-|U* = u*^ + (92 + § 4 ) 

= u* + (g-| + § 3 ) 
un+1 = ^n + At (A] + A2 )u^"^^/^ 


(XT ,X2> 

e Di 

(45a) 

(XT ,X2) 

€ §2 

(45b) 

(xi ,X2> 

e D 

(45c) 


No new boundary residual terms have been introduced for the explicit corrector 
(eq. (45c)), since this step is consistent with equation (1). Formation of the 
whole-step scheme requires care, since the functions g^ are discontinuous. 
Formulas derived for the contribution to the residual on D from the terms 
g^ in which operator commutativity was assumed did not correctly predict the 
computed results. For this reason, no operator commutativity was assumed on 
g^ in the following derivation. Therefore, formation of the scheme proceeds 
as follows. Multiplying equation (45b) by A-| and substituting equation (45a) 
for A-j u* yields 


A -| A2U >^+'' /2 = + (92 + 94 ) + A] (gi + 93 ) 


This equation may be solved by successively inverting the A-) and A 2 
operators: 

un+1/2 = a^^Ai'^uI' + A2^a7^ [(92 + 94 ) + Ai (gi + 93^ 

Substituting this into equation (45c) for multiplying the resulting 

equation by A-| A 2 produces 

A'|A 2 U*^'*"^ = A]A 2 U^ + At A*| A 2 (A"| + A 2 ) (A*|A 2 )“^u^ 

+ At A‘|A2(A-| + A2) A2^ A*| ^|jg2 + 94) A*| (g-| + 93^ 
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Commuting the operators A-|A 2 and (A-j + A 2 ) that are operating only on u” 
and simplifying the resulting equation gives the whole-step form; 


un+1 _ un 


un+1 _ yn' 


I yn+1 + yn\ ^^2 

— - (A, . Aj,)l I - — A,A2^^— 


AiA2(Ai + A2)A2^Ai' jj92 + 94> "*■ Ai (gi + 93)1 (x^ ,X2> ^ D (46) 


Figures 7 and 8 present the solutions obtained with LODBEC and LODPC, 
respectively. Note the minimuin residual of machine zero and the large maximum 
residual on It is interesting that the residual is machine zero on grid 

points immediately adjacent to ft for the LODBEC scheme for all time steps. 
(This is also true of the consistent LODFE scheme.) The boundary-error contri- 
bution to the residual moves in one extra row of grid points for the LODPC 
scheme, as shown in figure 8 - This is apparently the effect of the complicated 
operator on g in the last term of equation (46) . 

It is also interesting to compare the maximum residual from the ccxnputation 
with the predictions of equations (35) and (36) for the LODBEC scheme. For this 
scheme the steady-state residual is, fron equation (44), 


R 


(A-| + A2)u + 


A] (9i + 93) 


At 


(47) 


where the term proportional to u^ has been dropped, since it vanishes as 
steady state is approached. Substitution of equations (3) into equations (35) 
and (36) and substitution of this result into equation (47) gives, for the 
residual, 


^ At 

R = (A-I + A 2 )u A2fi (X-|,X 2 ) ^ w-| (48) 

h2 


R = R + — sin ttx 2 (x-j 4 X 2 ) ^ w-| (49) 

h2 


Thus, if R = 0, i.e., if the scheme has converged to steady state, then 


R = (A-j + A2 )u 


9 

TT-^ Sin TTxo 

h2 


(50) 
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and, hence 


R 


max 



(51) 


The computed results are compared with results predicted by equation (51 ) in 
table II for various values of x for a grid with h = 1/8. The computed and 
predicted results agree to within the spatial truncation error. The third 
column in table II gives the values predicted by equation (51) with A 2 f] 
evaluated discretely. The perfect agreement of the predicted and computed 
results confirms the validity of equations (33) to (36). 


REMDVAL OF SPATIAL INCONSISTENCY 

The first method for the removal of the spatial inconsistency investigated 
in the present work is the method of undetermined functions (MUD) as developed 
in reference 1 . The method of undetermined functions can be cast in two forms: 
one in which the boundary data are adjusted (MUDB) and another in which a cor- 
rection function is applied over the entire domain D (MUDF) . Understanding 
of MUD is essential to understanding the development of the general consistent 
scheme given later. 


Boundary-Point Method (MUDB) 

The basis for the MUDB LODBEC scheme lies in the equations defining g^^, 
i.e., equations (33) to (36), and the observation that since each of the steps 
in the method given by equations (42) is not consistent with equation (1), 
then u* is not necessarily an approximation to u at any time level. We are 
thus free to view f* as unknown boundary data which can be adjusted to yield 
a consistent whole-step scheme. In particular, if we set g-| = g 3 = 0, then 

* * 

equations (35) and (36) an be solved for f-] and f 3 : 


* 



-At 




n+1 


(XT ,X2) e 


(52a) 


f3 = f3 - At A2f3 = A2f3 (xi,X 2 ) e Y 3 


(52b) 


These modified boundary data are now applied on the first sweep for the calcu- 
lation of u*. Similar formulas for the calculation of f* for the consistent 
LOD forward Euler scheme (LODFEC) are presented in appendix A. 
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Field Correction Method (MUDF) 

A more general method of undetermined functions (MUDF) may be used to cor- 
rect for the boundary error in a way that allows use of the given boundary data 
in the calculations. For this method the LODBEC equations (eqs. (42)) are 
rewritten as 


A-j u* = u^ + At^A-|A 2 U^ 

(XT ,X2) e Dt 

(53a) 

A 2 U^"^^ = u* + q 

(X-] ,X2) e D2 

(53b) 


The unknown function q is introduced in equation (53b) for the purpose of 
absorbing the error in D due to the application of the given boundary data. 
The whole-step scheme becomes 


u^+l - un . 1 o ( 

- = (A-j + A2 )u^'^* - At^A-|A2 + A-|q (x-| ,X2) ^ D 


At 


At 


(54) 


Obviously in equation (54) we want A-j q = 0 for (x-) ,X 2 ) ^ D, We set 




q = 91 

(XT ,X2) e Yl 





) 

(55) 

q = 93 

(XT ,X2) e Y3 






The desired correction function q 

is obtained along each 

X 2 = Constant grid 

line (i.e., each Hj 

in D-j ) by solving 


Aiq = 0 

(XT ,X2) s Dt 


(56) 


along each such line subject to the boundary conditions given by equations (55) . 

Results of computations using the LODBEC scheme with the spatial inconsis- 
tency corrected by MUDB and MUDF are presented in figure 9. (Results of both 
methods are identical.) Note on the figure 9 residual map the random distribu- 
tion of machine zero residual over the entire grid and on the error map the 
maximum error at the middle of the grid. The results of these computations 
agree exactly with the results in table I. This agreement verifies that the 
scheme reproduces the exact solution to the difference equations generated with 
the Douglas-Gunn scheme for all three values of t. We thus term this scheme, 
i.e., LODBEC with MUD, complete consistent. 
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A SELF-CANCELING BOUNDARY-ERROR SCHEME 


Attempts to use a rigorous development of MUD for the LODPC scheme in the 
present investigation met with failure. The authors were able to develop a 
successful formulation of MUDF only by trial and error. It is thought that the 
reason for the failure to rigorously derive the proper field correction func- 
tion was the requirement to assume operator commutativity in order to obtain 
analytic formulas. Rather than presenting a detailed discussion of this situa- 
tion, we take here a more fruitful approach. Elimination of u* between 
equations (45a) and (45b) gives 


u'^ = AiA2u’’'''V 2 _ ^g2 + g2^ - A] ^gj' + 

Similarly, 

^n+l _ A-,A2nn+3/2 _ ^§ 2 ''’'’ + - A] + § 3 ''’^^ 

Substitution of equations (57) and (58) into equation (45c) gives 


(57) 


(58) 


^n^■3/2 _ 
At 


(un+3/2 ^ ijn+1 /2^ 

(Al + A2) 


At^ 

— At A 2 

4 


+ 





(59) 

Equation (59) indicates that the LODPC scheme is canplete consistent for the 
(n+1 /2) level data. In addition, note that no assumptions of operator 
commutativity were necessary to obtain equation (59). Ccxnputation with this 
method on the test problem produces results identical with results from the 
Douglas-Gunn method and from LODBEC with MUD. Note, however, that no modifica- 
tion of the boundary data had to be made to achieve these results. Also note 
that the n-level data converge to a different steady state than the (n+1 /2) 
level data. (See fig. 8 for the n-level solution to the test problem.) Thus 
the n-level data may be viewed as data that have been preconditioned by 
including the negative of the values of both the temporal and spatial inconsis- 
tency errors. Thus, the LODPC scheme, with the (n+1 /2) level data monitored 
as the solution, is a self-correcting LOD method which is complete consistent. 
To the authors* knowledge this result has not been observed before. 


26 



A GENERAL COMPLETE CONSISTENT LOD SCHEME 


The field method of undetermined functions (MUDF) may be used to develop 
a class of complete consistent LOD schemes. We will develop the algorithm 
first in two dimensions, demonstrate the accuracy with computational results, 
and then extend the method to an arbitrary number of split steps. 

We begin the development by writing the form of the desired whole-step 
formula 


r n - u^) 

= u^ + At(A*| + A2) + (1 - a)u^J - At^A-|A2 


(60) 


where 0 S a ^ 1 . The LOD scheme is then 


Al 


□* = Qn + ^ 


n+1 

<32 <34 j 


(xi ,X2) e Di 


(61a) 


^uH+1 = u* + + 93 '*’^ j (xi,X2) ^ D 2 


(61b) 


where 


Aj^ = (1 - a At A£) £ = 1 ,2 


(62) 


and with £ = are error terms due to the application of boundary 

data and are analogous to the g^^ of LODBE. Here is introduced to absorb 

the temporal and spatial error. Elimination of u* from equations (61) gives 


u^+l = qh + a At(A'| + A2 )u^"'‘^ - At^A-j A2U^'*'^ + + 94 "*"^^ 

+ A-j^g^'*'^ + g3^^^ (x-|,X2) G D u r (63) 

Now equation (63) is subtracted from equation (60) and the result is solved 
for Q^: 


QH 


u^ + (1 - a) At{A*| + A2 >u*^ + At^A-|A2U^ 

I ^n+1 , -^n+l\ ^ I -^n+1 ^n+l\ \ d 1 1 p 

- (92 + 94 1 - All 91 + 93 ) (X] ,X2) e d u F 


( 64 ) 
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By making use of equation (63) rewritten for this result can be stated as 


gn = gn -1 + At (At + A 2 )u*^ - ^ 92 ''’^ - § 2 ] - ^ 94 "^^ - 94 


- At 


''n+1 [ /^n+^ -vii 

91 -91+ 93 - 93 


(x-i ,X2) e D u r 


(65) 


For steady boundary data, equation (65) reduces to 


gn = gn 1 + At (A-| + A2)u^ (x-|,X 2 ) ^ D 


( 66 ) 


The scheme developed above is termed the LODQ scheme. For steady boundary 
data, the heretofore unspecified g functions in equation (65) may be dropped 
without introducing any error. For unsteady boundary data, however, neglect 
of the g function terms appearing in equation (65) produces an 0 (At) trun- 
cation error, as is now demonstrated- The terms due to g 2 and §4 need not 
be considered, since they will not affect the solution on the interior of the 
domain. A derivation completely analogous to that yielding equations (40) 
gives for g-| and 


Al 







(XT ,X2) e w-] 


(67a) 


Al 







(x^ ,X2) e W3 


(67b) 


Also, in the same way equations (35) 
for g-} and g 3 are obtained: 


and (36) were developed, expressions 


^ n+7 _n+l 
91 = fl 


f"! - 01 At A2f^ 


(XT ,X2) s Yi 


( 68 a) 


93^^ = - £3 - a At A 2 f 3 ^^ (XT,X 2 ) e Y3 ( 68 b) 

^ n“^~l ^ n-i”l 

For convenience, we now take f-j = f*| and f 3 = f 3 . Then equations ( 68 ) 
become 
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(XT ,X2) € Yl 


(69a) 




^n+1 . _n+l 

91 = -a At A2fi 


~ A ^n+1 

93 = At A2f3 


(xi »X2> e Y3 


(69b) 


Substitution of equations (69) into equations (67) then gives 


- g") = A2^fi''^^ - f") (xi,X2) e Yl 


(70a) 


Al 


^n+1 -^nN 
93 - 93j 


a At 

IJ' 


a At A 2 ( f 3 ^^ 



(x-| ,X 2 > s Y3 (70b) 


The reconstructed whole-step formula thus becomes/ with the use of equation (66) 
for Q, 


,n+1 _ „n 


At 


un f- -I (u'^+1 - u'^) 

— = (A-j + A 2 ) [au’^'*’^ + (1 ~ ot)u’^j - 0t2 At^AiA2 ^ ^ 


+ a 





(XT ,X2) € D (71 ) 


Now proceeding to the limit of vanishing mesh according to the law 
At/h^ = Constant shows that the scheme is consistent with equation (1 ) 
to 0(At). 

The canplete consistent LODQ scheme may be summarized as 


At u* = Q’’ 

(XT fX2) 

G 

Dl 

(72a) 

A2U^‘*'^ = U* 

(XT /X2) 

G 

D2 

(72b) 

gn+1 = gn + At (At + A 2 )u'^+'' 

(XT fX2) 

G 

D 

(72c) 
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where may be taken as u^. It should be borne in mind that for problems 

with unsteady boundary data# the scheme is first order in time regardless of 
the value of a. In order to achieve the possibility of second-order accuracy 
in time/ i.e./ to recover the whole-step formula given by equation (60) rather 
than that given by equation (71)/ one must evaluate the g-functions. Further/ 
it should be noted that operator commutativity does not have to be used to 
derive the equation for Q nor to develop the whole-step formula- 

The scheme given in equations (72) was run on TP2 with the results sum- 
marized in figure 10. The numerical solutions agreed to machine accuracy for 
all values of a and At with the Douglas-Gunn scheme results. 

The scheme can be readily extended to the case of m split steps as 


A”! u 

A2u(2) = uO) 


) (73) 

• » 

^u^+1 = u ) 

gn+1 = qh + At(A-| + A 2 + • . . + 


The scheme can be verified to be complete consistent by elimination of all 
intermediate values of u^^^ and Q. (The whole-step scheme for the 3-D case 
is given in appendix A.) The scheme is also first order in time for unsteady 
data regardless of the value of a. Again, operator commutativity does not 
have to be assumed in order to derive these results. 

A comment should be made about the stability of the LODQ scheme. As shown 
in appendix B, the method is unconditionally stable in two dimensions for both 
linear parabolic (heat conduction) and hyperbolic (advection) model equations 
for second-order central differences. In three dimensions/ however, the method 
retains unconditional stability for the parabolic case only. In the three- 
dimensional hyperbolic case, the method is unconditionally unstable for second- 
order central differences. This unfortunate result follows from the desired 
whole-step formula, which contains the destabilizing term ol^ At^6^u/6x] < 5 x 2 <Sx 36 t . 
One can eliminate this term from the whole-step formula, but to do so would 
generate a form of Q which would destroy the tridiagonal structure of one of 
the sweeps. This instability thus appears to be inherent in schemes for three 
dimensions which involve three spatially split implicit steps and two time 
levels, including the two-level Douglas-Gunn type schemes. 
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CONCLUDING REMARKS 




It is obvious fran the results presented in this paper that one must be 
very careful when applying locally-one-dimensional (LOD) type splitting schemes* 
Two types of errors present in such schemes have been explored in detail: the 

temporal inconsistency due to splitting and the spatial inconsistency due to 
application of Dirichlet boundary data. The temporal error has been shown to 
be of the order of the temporal truncation error of the scheme and, hence, does 
not have a large effect on the solution for explicit methods, in which one is 
forced to take small time steps from stability considerations. For implicit 
methods, however, the errors introduced by this inconsistency can be substantial 
when large time steps are taken. Care must therefore be taken to remove this 
inconsistency .in implicit schemes. 

The spatial inconsistency due to the application of Dirichlet data is some- 
what more troubling. As has been shown, this inconsistency produces 0(1) 
truncation error terms on grid lines adjacent to the boundaries and, hence, can 
have a substantial effect on solution accuracy even at small time steps At. 

For the test problem examined in the present work, the errors in the solution 
caused by this inconsistency proved to be acceptable (i.e., 0(h2), where h is 

spatial step size) at time steps of the order of the explicit stability limit, 
but such may not be the case for other problems. 

As a result of these observations a complete consistent LOD scheme has been 
developed. For steady Dirichlet boundary data, the method eliminates both iden- 
tified inconsistencies exactly; but, for unsteady Dirichlet data, the spatial 
inconsistency is made 0(At) rather than 0(1). This complete consistent 
scheme was developed in such a way as to produce a whole-step formula identical 
with that which is associated with the spatially split Douglas-Gunn method. 

As a consequence of this development the complete consistent locally-one- 
dimensional scheme (LODQ) can t>e considered an alternate interpretation of a 
general class of complete consistent ADI schemes which also includes the spa- 
tially split Douglas-Gunn method. All the strengths and weaknesses of this 
class of schemes are thus shared by the LODQ scheme. 

One of the principal shortcomings of the whole-step formula associated 
with the LODQ scheme is its unconditional instability for the 3-D linear 
advection equation. With the introduction of the method for development of the 
LODQ scheme, one has substantial freedom, however, in the choice of resulting 
whole-step formula. It may be possible to devise a whole-step formula which 
can be spatially split, complete consistent, and second order in space and still 
be stable for the three-dimensional advection equation (e.g., by introduction 
of another time level). These possibilities were not examined in the current 
investigation. 

One further comment should be made regarding the errors associated with 
LOD schemes. The time-split MacCormack-Baldwin method for the linear problem 
examined in this work reduces to the LOD Lax-Wendroff scheme given in appendix A. 
As shown in equations (A6) through (A7) this method contains the spatial incon- 
sistency and, as demonstrated in equation (A8) , also contains a temporal incon- 
sistency. Computational experience of the authors verified these formulas. 
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The inclusion of the cross-derivative terms in each split step of the MacCormack- 
Baldwin method for the full Navier-Stokes equations may circumvent the spatial 
inconsistency encountered in the present work. On the basis of the present 
results, however, one would anticipate the spatial inconsistency to appear in 
the implementation of the optimal splitting developed in NASA CR-1 59386. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
February 13, 1981 
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APPENDIX A 


OTHER LOP SCHEMES 

CONSISTENT LOD FORWARD EULER SCHEME 

In order to render the LOD forward Euler method complete consistent, a 
correction term similar to that in equation (42a) is introduced. Only the 
boundary method of undetermined functions (MUDB) for rendering the scheme 


complete consistent will be given. The algorithm 

is 


u* = Btu*^ + ^d 2 + d 4 ^ 

(XT ,X2) e D u r 

(A1 ) 

= B 2 U* - At^ATA2U^ + |d^ + d^ 

(XT ,X2) e D u r 

(A2) 

where A*|A 2 U^ = 0 for (xi,X 2 ) ^ and 




d" = - f* - At A2f * 

= f 2 - - At Aif^ 

d^ = - f 3 - At A 2 f 3 

d4 = £4 - £4 - At A-,£4 


(x-|,X 2 ) e Yl (A3a) 

(x-|,X2) e Y2 (A3b) 

(XT,X2) s Y3 (A3c) 

(XT ,X2) e Y4 (A3d) 


The whole-step £ormula is 


,n+l _ „n 


At 


— = (At + A2)u^ + At ^B2^2 (^1 »^2) ^ ^ 


(A4) 


In order to render the scheme complete consistent, we set d 2 = d^ = 0 by 
* * 

adjusting f 2 and according to the formulas 


APPENDIX A 


f2 = f” + ^1^2 (X1/X2) s Y2 


(A5a) 


£4 = £4 + At A-|£4 (Xi,X2) ^ Y4 (A 5b) 

Computational experiments conducted by the authors demonstrated the complete 
consistency of the above outlined method. 


LOD LAX-WENDROFF SCHEME 

An LOD scheme based on the classic one-step Lax-Wendro££ scheme can be 
written as 


u* = B] u*^ + - At^A-| A-| u*^ + ( e2 + ©4 


(X] ,X2) e D] 


(A6a) 


yn+1 = B2U* + - At^A2A2U* + ^e^ + e3^ (x-| ,X2> ^ D2 


(A6b) 


where 


e” = £^^^ - £l - At A2 £i - - At2A2A2£-| (x-| ,X2) ^ Yi 


(A7a) 


©2 = f2 - ^2 - At 1 At2AiAi£5 


(XI, X2) e y2 


(A7b) 


e" = £3'*’^ - £3 - At A2f3 - - At2A2A2f3 (xt ,X2) ^ Y3 


(A7c) 


©4 = £4 - £4 - At A-| f 4 - - At2A-|Ai£4 (x-| ,X2) ^ Y4 


(A7d) 
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The whole- step scheme is 


un +1 _ yjp 

At 


(Ai 


+ A2)u^ + - 


At(A*| + A2 )^u^ 


+ 


1 

2 


At2A-| A 2 (Ai 


+ A2)u*^ + 



A 2 A 2 A 1 A-| 


+ 


At"*^ ^82 + - At2A2A2^ ^02 + €4^ 


(xi ,X2) ^ D (A8) 


This scheme thus has both splitting and boundary errors. Computational experi- 
ence with this scheme reveals a steady-state solution which depends upon At^, 
even with the spatial inconsistency removed with MUD. 


LOD HYBRID SCHEME 

Another complete consistent LOD scheme in two dimensions is the hybrid 
scheme 


u* = 8*1 u^ + ^d 2 + < 34 ^ 

A2U*'+'' = u* + + 93^ 

The whole-step scheme is 
un+1 _ un 

^ = (Ai + A2 >u*^ + At 


(x-i ,X2) e Dt 
(XT ,X 2 ) e D2 


(un+1 _ un) 

A2 ;; (X] ,X2> e D 

At 


(A9) 


(AlO) 


(All) 


The scheme is thus observed to have no spatial or temporal inconsistency on the 
interior of the domain. Computational results confirmed this prediction. The 
scheme cannot be extended to three dimensions as a complete consistent scheme 
without adding terms and using MUD. 
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3-D LODQ SCHEME 

In three dimensions, the LODQ scheme has a split form very similar to the 
2-D case. The split steps are 


A-] u n ) = qU 

(XT ,X2,X3) 

e Di 

(AT 2a) 

A2u(2) = uO) 

(x-| ,X 2 ,X 3 ) 

e D 2 

(A1 2b) 

A 3 un +1 = u( 2 ) 

(xi ,X2,X3) 

G D 3 

(A1 2c) 

Qn+1 = Qn + At(Ai + A 2 + A 3 ) un+1 

(x-| ,X 2 ,X 3 ) 

G D 

(Al 2d) 


^ _ _ 

where A 3 and D 3 are defined analogously to Ajj, and Djj,, with 36 = 1,2. The 
whole-step scheme obtained by eliminating and Q is, for 

steady boundary data. 


un+1 _ 

At 


(Al 


+ A 2 



+ (1 



- At2(A-,A2 + A1A3 + A2A3) 


(un +1 _ 
_ 


+ At^A-| A 2 A 3 


(un+1 _ un) 

At 


(AT 3) 
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STABILITY OF LODQ SCHEME 
TWO DIMENSIONS 

In this section/ the Von Neumann damping ratios for the LODQ scheme are 
presented for the two-dimensional model equations 


9u 

9^u 9^u 


-f -■ - ■■ 

9t 

9x^ 9x^ 


(B1 ) 


8u 9u 9u 

9t 9x-| 9x2 


(B 2 ) 


For the two-dimensional parabolic equation (eq. (B1 ) ) , the damping ratio is 


1 + {a - 1 ) ( 3 i +62) + ot23iB2 

G “ ■ ■■ 

1 + ot(3] + &2) + 32 


(B3) 


For the two-dimensional hyperbolic equation (eq. (B2) ) , the damping ratio is, 
for central-space and backward Euler time differences/ 


G = 


n - a2mri2) “ (a - 1 ) i (m + TI 2 ) 
n - a2n-|n2) - oii(ni + H2) 


(B4) 


where i = and 

these formulas, \ g \ 


rijj, = At h ^ sin It is easily verified that in both 


^ 1 for all At provided that - ^ a ^ 1 . 


THREE DIMENSIONS 

In this section, the Von Neumann damping ratios for the LODQ scheme are 
presented for the three-dimensional model equations 
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3u 

3t 


9^u 3^u 3 2^ 

3 x^ 3 x| 9 x 2 


(B 5 ) 


3u 3u 3u 3u 

3t 3x-] 3x2 ^^3 


(B 6 ) 


For the three-dimensional parabolic equation (eq. (B 5 ) ) , the damping ratio is 

1 + (a - 1 ) (3i + ^2 ^ 3 ) + ^1^3 ^ 2 ^ 3 ) ot^3i3233 

G = — (B 7 ) 

1 + 01(3*1 + 32 + 63) + 0t2(3-[32 + ^1^3 ^2^ 3^ ^^^1^2^ 3 

1 

Again, from this formula, |g| = 1 for all At for - ^ a ^ 1. The damping 

2 

ratio for the three-dimensional hyperbolic equation (eq. (B6) ) is 

[T - a2(r)^ri2 + runs + n 2 n 3 r] + i |ja - 1 ) (ni +112+^3) - a^niri2n^ 

- a2(mri2 + niri 3 + n 2 H 3 ^ + i[ot(m + TI2 + TI3) - a^nin2n^ 

From equation (B8) we see that the scheme is stable provided 

|ja - 1) (m + ri2 + H3) - a3mri2n3j^ ^ [a(m + ri2 + ^ 3 ) ~ a^^in2n^^ 

This inequality can never be satisfied for At > 0 for any 01 in the range 
0 S a ^ 1 . Thus, the scheme is unconditionally unstable for equation (B6) . 


(B 8 ) 


(B 9 ) 
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TABLE I.- ERRORS IN SOLUTION TO TP1 AND TP2 


USING THE DOUGLAS-GUNN METHOD 

[t = 1 , 10, 1 ooo] 


h 

— 
1 ^Imax 

1 ^1 min 

1 ^ i max 2 

1/4 

2.8239880 x 1 0~^ 

1 .6879493 x 1 o'^ 

0.9037 

1/8 

7.293631 7 X i 0~3 

1 .601 0933 X 1 0-3 

.9336 

1/16 

1 .8392306 X l 0~3 

1 .211 291 3 X 1 0-3 

.9417 


TABLE II.- CALCULATED AND PREDICTED MAXIMUM RESIDUALS FOR 
TP 2 USING LODBEC METHOD 

[h - 1/s] 



1 ^1 max 

T 

Calculated by 
LODBEC 

Predicted by eq. (51) 
(analytic) 

Predicted by eq. (51) 
(discrete) 

1 

2.4358 

2.4674 

2.4358 

1 0 

243.58 

246.74 

243.58 

1 000 

2435.8 

2467.4 

2435.8 
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42 


Boundary conditions 


Analytic solution 


Test problem 1 (TPl) 



Boundary conditions 


Analytic solution 


Test problem 2 (TP2) 


Figure 2.- Boundary conditions and analytic solutions for the 
two test problems considered. 







Figure 3.- Numerical results - LOD forward Euler scheme applied to TPl . 



Numerical solution Error Residual 

(a) T = 1 , h = 0.1 . 



Numerical solution Error Residual 




U1 


(b) T = 1 0, h = 0 .1 . 

Figure 4.- Numerical results - LOD backward Euler scheme applied to TP2 



• Boundary points 
O Unknowns (treated implicity) 

® Other field points (treated explicitly) 



(a) Douglas-Gunn method: (b) lOD methods: X] line 

x-j line algorithm. algorithm. 


Figure 5.- Comparison of field points used in the line 
algorithms for the Douglas-Gunn and LOD methods. 
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Figure 6.- Propagation of the boundary residual onto the set 
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Numerical solution Error Residual 

Figure 9.- Numerical results - consistent LOD backward Euler scheme 
with MUDF applied to TP2. T = 10; h = 0.1 . 
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Figure 10.- Solution error of the LODQ scheme applied to TP2. 
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